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Abstract: The finite-difference time-domain (FDTD) method [1] is one of the most widely used full-wave 
electromagnetic simulation tools for solving Maxwell's curl equations due to its simplicity, 
straightforwardness and easy implementation. However, the Courant-Friedrich-Levy (CFL) stability 
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condition, i.e., Atrprp < 1/ (v 
limited the intense utilization of the FDTD. Thus, the development of accurate and efficient numerical 
algorithms for solving Maxwell’s equations 1s still an attractive topic in our computational electromagnetics 


(CEM) society. 


The precise integration time domain (PITD) method [2] is proposed for solving Maxwell’s equations, 
which breaks through the CFL limit and improves the computational efficiency. The fundamental idea of 
the PITD method consists of reducing Maxwell’s curl equations to a set of ordinary differential equations 
(ODEs) by approximating the spatial derivative with difference only, and solving the ODEs by using the 
precise integration (PI) technique with an accuracy of machinery precision. Compared to the 
convectional FDTD method, the PITD method posseses some striking features, which are summerized 
hereinafter. 

First, different with that of the conventional FDTD method, the stability condition of the PITD 


method can be expressed as Atpirp € V2l/ (v I F no? m Z) , where 1 =2" and N is a preselected 


integer. It can be clearly seen that the upper limit of the time-step size of the PITD method is determined by 
the preselected integer N and the spatial mesh size. Hence, the requirement of the time-step size can be 
satisfied by the alternative choices of the integer N . Meanwhile, it should be emphasized that the time-step 
size of the PITD method is much larger than that of the FDTD method. 

Second, although the numerical dispersion error of the PITD method is slightly larger than that of the 
FDTD method, the numerical dispersion performance of the PITD method is much better than other 
unconditionally stable algorithms, such as alternating direction implicit FDTD (ADI-FDTD) method. 
Furthermore, the numerical dispersion errors can be made nearly independent of the time-step size [3, 4]. 

Finally, the PITD method has advantages over the conventional Yee's FDTD method in using a large 
time step, and over the ADI-FDTD method in having high computational accuracy, respectively, and 
numerical dispersion errors of the PITD scheme can be made nearly independent of the time-step size. 
Therefore, we can use large time-step size to solve Maxwell's equations without increasing the 
computational error. 

Due to its excellent performance on computational electromagnetics, the PITD method attracts more 
and more attentions [2-7]. Recently, various low-dispersion and low-memory-requirement PITD methods 
have been proposed, such as the fourth-order accurate PITD [PITD (4)] method [3], the wavelet Galerkin 


scheme-based PITD (WG-PITD) method [4], the leapfrog scheme-based PITD (L-PITD) method [5], the 
unified split-step PITD (SS-PITD) method [6], the compact PITD (CPITD) method [7] and so on. With the 
rapid development of the computer hardware technology, the PITD method still has 
considerable development potential in the future. 
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CEM map “ 


Hybridisation to 
solve large and 
complex 
problems 


ELECTRICAL SIZE 





* COMPLEXITY OF MATERIALS 
source: http://feko.info/product-detail/numerical_methods 
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FDTD timeline 
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FDID major technical paths 


VY Absorbing boundary conditions 

e Mur; Engquist-Majda; Berenger PML, UPML, CPML 
Numerical dispersion 

e High-order space differences; MRTD; PSTD 
Numerical stability 

e ADI techniques; PITD; One-step Chebyshev method 
Conforming grids 

e Locally/globally conforming; Stable hybrid FETD/FDTD 
Digital signal processing 

e Near-to-far-field transformation 
Dispersive and nonlinear materials 

e Isotropic/anisotropic dispersions; Nonlinear dispersions 
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Multiphysics 
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PITD timeline 


PI techniqe for CEM (Ma) 2000 

















ybrid PITD-FDTD 
Krylov subspace-based PITD 


PITD monograph 
(Sci. Press, Ma) 
Unified split-step PITD (Ma) 


Compact PITD (Ma); 
Split-step-scheme-based PITD (Ma) 


Leap-frog scheme (Ma); 
Lossy media (Ma) 


High order PITD (Ma) 


PI technige for 
LTI ODES (Zhong) 
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Sub-domain PITD (Ma); 
a L 2010 "wavelet Galerkin PITD (Ma) 


[ETE Arad 





Gn Ab en ICEF2016 Xi'an -- PITD 6/35 


(C23 


ICEF 2016 Uvervie\ Methods EXampies 








Maxwell's equations 


And God said: 
CE — el . V x H, 
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OH 1 
d ER. E 
Ot VEE 





and then there was light. 


J. C. Maxwell 


“From a long view of the history of mankind the most significant event of 
the nineteenth century will be judged as Maxwell's discovery of the laws 


of electrodynamics." --- Richard P. Feynman 
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Discretization and Yee cells 





FDTD PITD 
e Finite difference in space e Finite difference in space 
e Finite difference in time e ODEs in time 
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Updating equations / ODEs 


FDTD PITD 
(R 4 F) Kort — (R — F) X" 4 fo+1 dX(t) ) = MX(t )+ F(t) 


dt 
oe) Gk - 3 


where 





å E” € M — matrix containing material 
Atc | H^*1/2 | properties and the discretization of the 
curl operators 
@ Doom — diagonal matrices 


containing €, fl, Te, Om for each cell 


€ f — sources 


@ K — arises from the discretization of 
the curl operators 


e f"*! — sources 
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Formal solution to & = MX + f(t) 


e Analytical form 


X(t) = exp(Mt)X(0) + [ exp[M(t — s)|f(s)ds 


e Recursive form 
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Ways to compute eM^t 


SIAM REVIEW 
Vol. 20, No. 4, October 1978 


@ Society for Industrial and Applied Mathematics 
0036-1445 /78/2004-0031501.00/0 


NINETEEN DUBIOUS WAYS TO COMPUTE 
THE EXPONENTIAL OF A MATRIX* 


CLEVE MOLER?+ AND CHARLES VAN LOAN? 


SIAM Review 


© 2003 Society for Industrial and Applied Math 
Vol. 45, No. l, pp. 3-49 


Nineteen Dubious Ways to 
Compute the Exponential of a 
Matrix, Twenty-Five Years Later* 





Cleve Moler* 
Charles Van Loan? 
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Series methods 

ODE methods 

Polynomial methods 

Matrix decomposition methods 
Splitting methods 


Krylov space methods 
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PI technique to compute T — eM 


© compute T = (eMry" - (14 T5)! -- to be contd. 
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PI technique to compute T — eM^* (contd.) 
© compute T = (emry" - 14T.) 
N—1 2N— 1 


T= (1+ Ta)? (4T. riri) a. 


with (I-- Tj)? 2 1- 2T, 4 T; x T, 


Pseudo code 





do í = 1,..., N 
T; 2T;+ T, x Ta 
end do The algorithm holds the 
T —I-T; computational precision. 
©MET 
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tn+1 
Evaluating T| e *Mf(s)ds & the recursive formula 
tn 


(1) Analytical solution under the linear representation of f(t) 


e 15 order Taylor approximation of f(t), i.e., 


f = fo + (t = tn)f1, te [Es tid) 
tn+1 
e integrate T| e *Mf(s)ds analytically 
tn 
e recursive form for X,11: 


Xn+1 =T [Xn + MT fo + Mtf )] 
—M [fo + MT + f At] 
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tn+1 
Evaluating T| e *Mf(s)ds & the recursive formula (contd.) 


tn 


(2) Recursive formula by using the Gaussian quadrature rule 


e three-point Gaussian quadrature 


e recursive form for X n+1: 
Xs 7 TX, i 3 etes) MAH ai + (1 i v35) At/2| 
3 (1-V35)Mat2, È È (1 + 3/5) At/2| 
18 
+ at (tn + At/2) 
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Remarks 


X 
Z = MX +f, X = [E,H]' 


@ finite difference in space, but differential in time 
@ scaling and squaring for exp(MAt) 


© T,<2T,+ T; x T, to guarantee the computational 
precision 


© Gaussian quadrature for the excitation term 


6 non-staggered E and H in time 


Precise Integration 
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Numerical stability condition 


1 
FDTD CFL criteria: Atupper = 


UNO eal ge + ay? + Az? 
PITD 


Stability conditions vary for different orders of Taylor approximation: 
e 18¢/2"¢_order: unstable: 
o 3-order: At < V3/20JAtERIP — V3/2 aig 


© Athorder: At < /2¢AthOLD — 2H A HEP Almost unconditionally 
oun 


LPer 





stable for large N. 
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Numerical dispersion analysis 
FDTD numerical dispersion relation: 
We/c* = Wo + W? + WÈ, 


— sin (ky y|z4x|y|z/2) W, = sin(wAt/2) | 


where W.y]z = Ax|y|z/2 At/2 


|ylz 


PITD 


2 
— EG) _ (App = Nérrp/3!) 
È 1 + Abrrp/2! — Abrrp/4! 


where Apitp = cat [We T W? T WÈ. 
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Numerical dispersion analysis (contd.) 
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Source and boundary conditions 
Source conditions Boundary conditions 


e Hard sources 


e Plane waves & TS/SF technique 


e Engquist-Majda ABC 
e PMLs 
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Remarks 


Characteristics of the PIT D method: 


"4 
"4 


preselected N determines the upper bound of At"!"D 


PITD FDTD 
bound boun 


slight worse numerical dispersion compared with that of the 
FDTD method 


numerical dispersion can be independent of At 


technique paths of the FDTD method can be learned 


“Stones from other hills may serve to polish the jade of this one.” 
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Improved methods 


e Fourth-order PITD [PITD(4)] method 

e Wavelet Galerkin PITD (WG-PITD) method 
e Leapfrog PITD (L-PITD) method 

e Compact PITD (CPITD) method 

e Hybrid PITD-FDTD method e 


e Krylov subspace method O EN 
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Improved methods — PITD(4) ! 


4*^- order spatial difference scheme is used as: 








ôu; 1 f1 27 2 
Ox  Ax|24 (ui-3yo — Uj+3/2) 94 (ui-iyo — U;+1/2) +0 | (Ax) 
E E 1.000 
3 £ oss| | Numerical dispersion 
E —+— FDTD 
2 > PITD(O) S 0.990 
— PITD(4) 
—-— FDTD(2,4) 
" " at "o ka ^" m 001 0.1 1 10 100 1000 10000 


Angle of propagation, deg Posi abes 


‘EEE T-AP, 59(4), 2011: 1311-1320. 
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Improved methods — WG-PITD method ? 


Discretization form of Maxwell equation(s) in space: 





d Exlj41/2,m,n oli+1/2,m,n Hzli41/2,m+i+1/2,n Hylja1/2,m,n+i+1/2 
—— n S oL Edrb2,ma * 218i rata: a 
dt Elj41/2,m,n ; #1141/2,m,n AY El14+1/2,m,n AZ 
ao døy(X) 
where aj = | , 4 0-;(x)dx. 


Numerical dispersion S11-parameter 


IS, dB 
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Relative error of numerical phase velocity 
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Improved methods — Hybrid PITD-FDTD method 3 


How to handle multiscale problems with fine geometrical features? 


e Subgriding in FDTD — At depends on AXmin ) 
e PITD — need to compute eM^*. put At can be relaxed 
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Improved methods — Krylov space method ^ 
Recursive form of the PITD method: 


Must s eM ae > aje BAN (t) + Y At) 


I 


number of nonzero 
Denseness Memory cost (MB) 


elements 
M 1558 0.0039 0.003 
ee 370482 0.9334 0.76 


Evaluate eM^* explicitly? X —> Estimate e™4‘v directly. Y 
p y y 


“Please refer to OC2-6 (Tue. 08:30-10:00, Function Room 3, 3F) 
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Direct estimation of eM^ty 


€) m"-order Krylov subspace: K"(M, v) — span(v, Mv,..., M" 1y) 


@ Arnoldi process 
oe Vm= [V1,V2,--. Mal” — orthogonal basis of K” 


e H, « V/, MV, - matrix generated during the Arnoldi process 
© M^v ~ Vme” ^V Iv = Vme" ^ei, e —[1,0,0,...,0] e R"*! 









x 10” 
3f A —— Analytical solution | 4 
3 # Proposed approach 
al. F A FDTD | M 
TA . emory 
E | 1 | CPU time (s) cost (MB) 
S 1 Krylov-PITD 23.78 0.30 
å FDTD 137.28 1.34 


DE 








1 1 1 
5 10 15 
Time/ns 


“Please refer to OC2-6 (Tue. 08:30-10:00, Function Room 3, 3F) 
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Rectangular cavity 


fana. (GHz) 


3.125 
4.881 
5.340 
7.289 
7.529 


FDTD scheme 


At= 1ps 
f (GHz) rel. err. 
2.983 4.54% 
4.750 2.68% 
5.450 2.06% 
7.333 0.60% 
7.567 0.51% 


Examples 


ADI-FDTD scheme 


At = 60 ps 
f (GHz) rel. err. 
2.900 7.20% 
4.650 4.73% 
5.580 4.49% 
6.817 6.92% 
7.000 7.03% 


@ relative error increases as the time-step 


increases for the ADI-FDTD method 


@ relative error is independent of the time 


step for the PITD method 


relative error(%) 
co 


eo 
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Una point of FDTD 


- 


PITD scheme 

At = 60 ps 
f (GHz) rel. err. 
2.983 4.54% 
4.750 2.68% 
5.450 2.06% 
7.333 0.60% 
7.567 0.51% 





—&- PITD 
—+ ADI-FDTD 
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Microstrip low pass filter ? 






S 2.54 mm 
ID 


IS, sail dB 


—+—|S,,|, FDTD 
—e—|S,,|, FDTD 
—+— |S,,|,L-PITD 
——|S,,|,L-PITD 





8 10 12 14 16 18 20 
Frequency (GHz) 


Comparison between the FDTD method and the L-PITD (Leapfrog PITD) method İ 


Ax Ay Az At memory CPU time 
FDTD 041mm 0.26mm 042mm 0441 ps 24.44 MB 1024 s 
L-PITD 0.41mm 0.26mm 042mm 0.884 ps 248.4 MB 851 s 
T Simulations were performed on Intel® Core™ Duo CPU T8100 2.10 GHz PC. 


IEEE MWCL, 22(6), 2012: 294 - 296. 
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Small antenna 


Coarse Grid 


Analytical solution 
(FDTD) T 


* Proposed approach 


å 








Subgrid 
(Krylov) 





Dipole 
Atenna 

















270 





Comparison between the FDTD method and the Krylov-PITD method 


CPU time (s) Memory cost (MB) 
Krylov-PITD 23.78 0.30 
FDTD 137.28 1.34 


Gn Ab n ICEF2016 Xi'an -- PITD 30/35 


(C23 


ICEF 2016 Outlooks 








Remarks 


W almost unconditionally stable 
V relative large time step size can be used 
PI technique maintains the computational precision 


relative error independent of the time-step 


NN S 


hybrid PITD-FDTD technique suitable for multiscale problems 


w memory cost can be relaxed by using the Kyrlov space method 


“Sey NFAT è” — GUI) 


“One flaw cannot obscure the splendor of the jade.” 


— Book of Rites w Meaning of Interchange of Missions twixt Different Courts 
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Outlook 


Future work : 
e Sub-domain technique 
e Parallel computing technique 


e Extend to complex materials 


Future prospects : 


e Nanophotonics and nanoplasmonics. Ultimately, combination 
of quantum and classical electrodynamics 


e Multiphysics 
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Thank you! 
Q & A. 


Recruitment Ad 


Our Group for Advanced Electrical Technologies (GAET) is 
seeking for highly self-motivated students with the solid scientific 
strength in one of the following areas: 


e Electromagnetics (theory and computation) 

e Metamaterials and plasmonics (graphene plasmonics) 
e Wireless power tansfer 

e Power electronics 


Please visit http: //tydong.gr.xjtu.edu.cn/ for details. 





